
# y = x^2
def f(x):
    return x*x

# trapzoid integration of f(x) on [a,b], with n divisions
def trapzoid1(a, b, f, n):
    s = 0
    h = (b-a)/n # 步长
    for i in range(n): #range(n) = [0,1,...,n-1]
        s += ( f(a + i*h) + f(a + (i+1)*h) ) * h / 2
    return s

# trapzoid integration of f(x) on [a,b], with error <= precision
def trapzoid2(a, b, f, precision):
    n = 10
    error1 = 1000 * precision
    Deltan = 1
    Delta2n = 1
    while( error1 > precision ):
        NIn = trapzoid1(a, b, f, n)
        NI2n = trapzoid1(a, b, f, 2*n)
        error1 = (NIn - NI2n)/3
        #print(" n = ", n, " integral = ", Delta2n)
        #print(" error1 = ", error1, "precision = ", precision)
        n = 2*n
    return NI2n, (NIn - NI2n)/3 

x = trapzoid2(0, 1, f, 0.00001)
print( "integral = ", x[0], " error = ", x[1] )
